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Abstract. We develop an analytical model for ultracold atom- ion collisions using the 
multichannel quantum-defect formalism. The model is based on the analytical solutions of the 
r -4 long-range potential and on the application of a frame transformation between asymptotic and 
molecular bases. This approach allows the description of the atom-ion interaction in the ultracold 
domain in terms of three parameters only: the singlet and triplet scattering lengths, assumed to 
be independent of the relative motion angular momentum, and the lead dispersion coefficient of 
the asymptotic potential. We also introduce corrections to the scattering lengths that improve the 
accuracy of our quantum-defect model for higher order partial waves, a particularly important result 
for an accurate description of shape and Fcshbach resonances at finite temperature. The theory 
is applied to the system composed of a 40 Ca + ion and a Na atom, and compared to numerical 
coupled-channel calculations carried out using ab initio potentials. For this particular system, we 
investigate the spectrum of bound states, the rate of charge-transfer processes, and the collision 
rates in the presence of magnetic Feshbach resonances at zero and finite temperature. 
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1. Introduction 

The successful realization of systems combining ultracold atoms and ions is stimulating an 
increasing interest in the physics of atom-ion collisions at very low collision energy. The first 
experiments were realized for clouds of atoms and ions stored in hybrid dual charged-neutral traps 
at niK temperatures [1,2]. More recent experiments performed with atoms in the micro and nano 
Kelvin temperature domain consist in immersing single ions in a Bose-Einstein condensate [3-5]. 
Understanding atom-ion collisions in the quantum regime is an essential elementary step before a 
more complete many-body description of these systems can be developed. While atom-ion collision 
properties are well known for high collision energies [6, 7], a theoretical description in the ultracold 
domain is still under development [8—11]. Atom-ion and neutral atoms scattering are significantly 
different, mainly due to the relatively long-range character of the polarization interaction between 
atom and ion. Beyond the fundamental interest of their collision properties, these systems are also 
very attractive for quantum information processing [12, 13]. Hybrid architectures may profit from 
advantages offered by both ionic and atomic species, namely a short computation time for charged 
particles and a long-coherence time for neutral atoms. Ultracold charged-neutral systems are also 
expected to exhibit interesting many-body effects, including for instance nontrivial modifications of 
the condensate wave function in the presence of ionic impurities [14] or the creation of mesoscopic 
size molecular ions [15]. 

In a previous paper [10] we have studied the basic properties of atom-ion scattering and 
bound states, using 40 Ca + — 23 Na [16] as a reference system. We have developed an effective 
atom- ion collision model by applying the multichannel quantum defect theory (MQDT) [17- 
19] to the polarization potential, which scales as r -4 at large atom-ion distance r. We 
have verified that MQDT model predicts very accurately all atom-ion scattering properties at 
ultracold temperatures by comparing to numerical close-coupled calculations performed on ah 
initio 40 Ca + — 23 Na potentials. In the literature, the MQDT approach has already been applied to 
describe scattering and bound states in electron- ion core [17], electron-atom [20] and neutral atom 
systems [21]. Our model of atom-ion collisions combines knowledge of the analytical solutions for 
the r -4 asymptotic potential [20,22,23] with the idea of a frame transformation [21,24-28]. The 
latter, applied at small distances allows a reduction of the number of quantum-defect parameters 
needed to represent the effect of the short-range interaction potential. 

In [10] we have applied MQDT to the study of magnetic Feshbach resonances and of radiative 
charge exchange processes, including the effects of Feshbach and shape resonances. In this 
paper we present a detailed derivation of our atom-ion model, based on the Mies formulation 
of MQDT [19]. We begin by introducing the quantum-defect approach, and discuss the properties 
of analytical solutions of the r -4 asymptotic potential. We derive asymptotic expansions of the 
model quantum-defect functions, valid for energies below a characteristic quantity E* defined in 
terms of the atom-ion interaction strength. For energies above E*, an efficient numerical algorithm 
to determine the characteristic exponent and the solutions of the Mathieu equation needed for the 
MQDT is proposed. Knowledge of the quantum defect parameters is crucial for understanding 
atom-ion collisions in the ultracold domain, and allows us to derive several non trivial analytical 
results concerning the behavior of weakly bound states, s-wave magnetic Feshbach resonances, 
and radiative charge-transfer probabilities. In comparison to the model presented in Ref. [10] 
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we introduce an angular-momentum- dependent correction that improves the accuracy of MQDT 
for high order partial waves. This modification allows the regime of applicability of our model 
to be extended up to mK temperatures. Using this revised model and numerical close-coupled 
calculations, the effect of finite energy Feshbach resonances on elastic and charge-transfer collision 
rates is investigated as a function of the magnetic field. 

In the present paper we also include details on the semiclassical description of radiative charge- 
transfer rates. A combination of the semiclassical approximation and MQDT scaling functions 
allows us to represent the charge-exchange loss rate as the product between the classical Langevin 
result and a universal quantum correction dependent only on a single parameter, the short-range 
phase. In [10] we discovered that the main product of a charge-transfer process for 40 Ca + — 23 Na are 
molecular ions. Here we present a more detailed analysis by analyzing the population of vibrational 
states formed after the electron transfer. A semiclassical formula describing this distribution is 
derived. The very good accuracy of our MQDT analytical results is tested by comparison with the 
numerical solution of the close-coupled Schrodinger equation for the 40 Ca + — 23 Na system. 

The paper is organized as follows. Section II introduces the MQDT formalism. Section III 
discusses the properties of scattering and bound states for a single channel potential, based on the 
analytical solutions for the polarization potential. In particular, section III. A derives the analytical 
solutions, section III.B presents the results of the scattering states and analyzes the applicability 
of the semiclassical description. Section III.C presents the analytical results for weakly bound 
states and section Section III.D studies the low-energy behavior of the MQDT functions. The 
idea of the frame transformation is discussed in Section IV, using a particular example of an alkali 
atom and an alkali-earth ion. Section IV.A generalizes the frame transformation to the presence of 
magnetic field, and Section IV. B derives the corrections to the quantum-defect matrix for higher 
partial waves. Section V describes the radiative charge transfer: the perturbative approach based 
on the Fermi golden rule, and the semiclassical description of the charge-transfer probability at 
short range. The theoretical model is applied to the 40 Ca + — 23 Na system in Section VI. Section VII 
presents the conclusions, and five appendices give some technical details related to the derivation 
of the analytical solutions and the multi-channel calculations of the radiative charge transfer. 

2. Multichannel formalism and quantum defect theory 

In this section we summarize the basic formalism of the multichannel quantum defect theory, 
adopting the formulation introduced by Mies in Ref. [19]. Atom-ion collisions are described by a 
iV-channel close-coupled radial Schrodinger equation 



Here fi = mim a / (mi + m a ) denotes the reduced mass, and W(r) and F(r) are N x N matrices 
representing the interaction and the radial solutions, respectively. A general solution to the N- 
channel scattering problem is given by a set of N linearly independent wave functions 




(1) 



N 




(2) 
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where |$,) are channels states describing the internal spin degrees of freedom, l^-(r) denotes the 
angular part of the solution (spherical harmonic) for the channel j. The interaction matrix is 
asymptotically diagonal 



where E°° are the threshold energies for the molecular dissociation, £i is the partial wave quantum 
number of channel i, C4 = ae 2 /2 with a denoting the static dipolar polarizability of the atom and 
e is the ion charge. Here we neglect the contribution of the higher-order dispersion terms to the 
long-range potential, which give relative small corrections in atom-ion scattering [10]. 

Given the total energy E, the channel states can be classified as open for E > E°° or closed 
for E < E?° . In the former case the asymptotic wave number k{ = ^Jl\i{E — Ef°)/h 2 is real and 
positive, while in the latter is purely imaginary ki = e m ' 2 \ki\. The solution matrix F(r) can be 
split into blocks 



corresponding to N open and N c closed channels. Imposing appropriate asymptotic boundary 
conditions on the closed channel components, Eij — y (r — > 00) for i = N + 1, . . . , N, the 
physically meaningful part of F(r) is contained in the block (N x N a ). The observable properties 
of the atom- ion system result from the asymptotic behavior of the open-open block F OQ (r), which 
at large distances yields the reactance matrix K OQ 



Here, Jij{r) — > 5ijSin(kir — ii / K/2)/y/ki and Nij(r) — > — cos^r — £j7r / '2) / \fk~i (r — > 00) 
exhibit asymptotic behavior associated with the spherical Bessel functions je(kr) and ne(kr). The 
constant matrix A DO depends on the boundary conditions at r — > 00. In particular, the choice 
A 00 = (1 — iK 00 ) _1 corresponds to the usual incoming-wave boundary conditions, where the 
amplitude of the outgoing wave is determined by the scattering matrix S DO = (l+iK 00 )(l— iKoo)" 1 . 

The basic idea of the quantum-defect theory is to introduce a set of parameters describing 
the short-range behavior of the wave function, that weakly depend on total energy E, and can 
be used to predict the system properties as E crosses the dissociation thresholds of the individual 
channels. As we will show later, in the ultracold domain the quantum defect parameters for atom- 
ion collisions are also weakly dependent on the relative orbital angular momentum I and can be 
taken as constant, at least for the lowest partial waves. 

We now specialize the formulation of MQDT developed by Mies for atomic collisions [19] to our 
atom-ion system. The starting point is the choice of a set of reference potentials {Vj(r)} that should 
reproduce the asymptotic behavior of the interaction matrix at large distances Vj(r) r — £ Wjj(r) 
but can otherwise be arbitrary. One associates to the reference potentials Vi(r) a pair of linearly 
independent solutions fi(r) and c/i(r) that have short-range WKB-like normalization 




(3) 




(4) 




(5) 



fi(r) = tti(r) sin/3i(r), 
9i{r) = oii(r) cos/3;(r). 



(6a) 
(6b) 
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The amplitude cti(r) fulfills the inhomogeneous Milne equation: [d 2 /dr 2 + ki(r) 2 ]ai(r) 



with the local wavevector ki(r) = y / 2fi(E — Vi(R)/h, while the phase df3i/dr = 1/af. Since 
the reference potentials {Vj(r)} reproduce the asymptotic behavior of W(r), the exact solution 
to Eq.(l) can be expressed at large distances in terms of pair of functions f(r) = {Sijfi(r)} and 

F(r) ™ [f (r) + g(r)Yl A (7) 

Here, Y is the quantum-defect matrix that play a central role in the MQDT analysis. In contrast 
to the S and K scattering matrices, Y remains analytic across the thresholds, and has in general 
only a weak dependence on energy. 

The observable properties depend on the asymptotic behavior of the solutions at large distance. 
It is therefore convenient to introduce another pair of solutions, fi(r) and gi(r), together with the 
physically well-behaved solution 4>i(r) for the closed channels, that have energy-like normalization 
as r — )• oo 

fi(r) = k; 1/2 sm(hr - £ i7 r/2 + E > E°° (8a) 

9i (r) 9* kr 1/2 cosfor - ii-K/2 + E > E™ (8b) 

Mr) = \k; 1/2 e~^ r , E < E°° . (8c) 

The factor | in the last equation is introduced for convenience in order to simplify the value 
of Wronskian [19]. Here, £j is the scattering phase shift for a channel i. The MQDT functions 
Ci(E), tanAj(E') and v^E) relate the solutions (8) to (6), and are defined as follows 

f t (r) = C7\E)Ur), E > E°° (9a) 

9i (r) = Q(E) [^(r) +tw\i{E)fi(rj\ , E > E°° (9b) 

Mr) = ftfi(E)^cosUi(E)fi(r) - sin Vi(E)gi(r)] , E < E°° (9c) 

Now, when all the channels are open: Vj-E > E°°, the solution F(r) at r — > oo can be expressed 
in terms of f(r) = {Sijfi(r)} and g(r) = {S^g^r)} 

F(r) ™ [f (r) + g(r)R(£)] A (10) 

Using the relationships (9a)- (9c) one can show that 

R(E) = C~\E) [Y- 1 (E)-tan\(E)]~ 1 C(E) (11) 

where C(E) = {SijCi(E)} and X(E) = {5ij\i(E)}. At high energies the WKB approximation is 
valid at all distances, and the functions /, g become identical to / and g, respectively. Therefore, 
in this situation the MQDT parameters behave like C{(E) — » 1 and tanAj(E') — > 0, and Eq. (11) 
reduces to R(E) = Y(E). 

Next, with the help of Eq. (5) one can relate the reactance matrix K and the scattering matrix 
S to the matrix R 

K = [sin(£) + cos(£)R] [cos(£) - sin(£)R]~ 1 , (12) 
S = e^[l + 2R][l-iR]" 1 e^, (13) 
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Table 1. Characteristic distance R* and characteristic energy E* for some selected atom-ion 
systems. 

R* (units of a ) E*/h (kHz) 

40 Ca+ + 23 Na 2081 28.56 

40 Ca+ + 87 Rb 3989 4.143 

135 Ba+ + 87 Rb 5544 1.111 

172 Yb+ + 87 Rb 5793 0.9313 



with = {SiME)}. 

The current derivation has to be modified in the presence of closed channels: E < E°° for 
i = N + 1, . . . , N. In this case, we impose the requirement that the wave function of the closed 
channels decays exponentially for large r, i.e. the closed channel wave functions are proportional 
to <fti(r). This results in the renormalization of the open-open block of the quantum-defect matrix 
[19] 

Y 00 = Y 00 - Y oc [tan(i/cc) + Y cc ] 1 Y co , (14) 

where u(E) = {6ijVi(E)} The scattering matrices can now be calculated from Eqs. (11)-(13), 
applied only to the open-open block, where one substitutes Y GO in place of Y OQ . Finally, when all 
the channels are closed (N c = N), the wave functions of all the channels must be proportional to 
4>i(r) at large r, and the energies of the bound states are determined by the condition 

| Y(E)+ tan i/(£)| = 0. (15) 

3. Long-range atom-ion interaction 

3.1. Analytical solutions 

Here we focus on the Schrodinger equation for a single channel where we include only the long-range 
part of the atom-ion interaction and the centrifugal barrier for partial wave £ 

In the following we will work in dimensionless units, where the length is expressed in units of 
R* ee v /2C 4 /i/n 2 and energy in units of E* = ti 2 / [2fi(R*) 2 }. Table 1 presents the characteristic 
lengths R* and energies E* for some sample combination of alkali atoms and alkali-earth ions. Fig. 1 
shows the long-range atom-ion potentials for the lowest partial waves, where the squares mark the 
top of centrifugal barriers occurring at r max = y/2/y/e(e + i)R* with energy E max = ^£ 2 (£ + l) 2 E* . 
The Eq.(16) in characteristic units of R* and E* takes the form 



This equation can be solved analytically [22,23,30] by substituting F{r) = ?/;(r)r 1//2 and 
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Figure 1. The long-range atom- ion potentials for the few lowest partial waves scaled by the 
characteristic distance R* and the characteristic energy E* . 



r = e z E 1 / 4 , which yields the Mathieu's equation of imaginary argument 



dz 2 



[a- 2gcosh2^]^ = 0. (18) 



with a = {I + |) 2 and q = \fE. Although the Mathieu's equation is well known in 
mathematical physics (see e.g. [31,32]), we include a brief derivation of its basic properties in 
the Appendix Appendix A. The Appendix also discusses the basic methods we applied for its 
numerical solution. 

We denote two linearly independent solutions to Eq. (17) by T v (r) and T_ u (r), where v is the 
complex number called the characteristic exponent (see Appendix Appendix A for definition). In 
view of the asymptotic properties of T„(z) discussed in Appendix Appendix A), we can construct 
two solutions f(r) and g{r) defined by Eqs. (6) 

/>) =AMT v {r) + A_MT-u{r) (19) 
g(r) =A u (<p + 7r/2)T„(r) + A- V (<j> + x/2)T- u (r), (20) 

where 

a <a>\ sin(0-z/7r/2 + 7r/4) 

AuW) = —, — v (21) 

sm(7rz/) 

and is some parameter that can be interpreted as the short-range phase. As can be easily verified, 
f(r) and g{r) exhibit at distances r C 1 semiclassical behavior given by Eqs. (6) with a(r) = r 
and 0{r) = -1/r 

f{r) ^ rsin(-l/r + 0), (22a) 
g(r) ^4 r cos(-l/r + 0), (22b) 



The asymptotic behavior given by Eqs. (22) can be also obtained by solving directly Eq. (17) 
with centrifugal and energy terms neglected. Since r sin(— 1/r + 0) is an exact solution of (16) for 
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£ = and E = 0, it can be used to derive a simple relation between the short-range phase and 
the s-wave scattering length 

a = R* cot 0, (23) 

that follows from the asymptotic behavior: F(r) — > r — a (r — > oo) at E — 0. 
At large distances, the solutions f(r) and g(r) behaves according to 

f(r) r ^C u (4>) y fqrj i {qr) - D u (<j))y/qrnt(qr), (24a) 

g(r) ™CV + |j y/qrj e (qr) - D v + |j y/qrn e (qr), (24b) 

where 

O(0) =4(^)171, cos 77 - (-l)M„ Jy (0)m_^sin?7, (25a) 

Dv(4>) —(—1)A- V ((j))m- V cos 77 — A v {i$)m v sin 77. (25b) 

Here, 77 = |(^ — £ — |), m„ = S v {A/q) v , and Sj, is a function of v, that is defined in terms of a 
continued fraction (see Appendix Appendix A for the definition). 

Making use of the asymptotic behavior (24) one can find the phase shift £ 

tan£ = A,(0)/O,(0), (26) 
and all MQDT functions defined in (9a)- (9c) 

C(E) =0,(0)/ cos £ (27) 

tanA(E) = C-^taiT^f-O (28) 

tan ^ ^ ^(0 +f)+ ^(0 + f)^(x/4)- (29) 



where x = V~E, and £ is the phase shift of a second short-range normalized solution g(r) 

tang = — , 30 

6^(0 + 7T/2) 

Finally we point out that formulas (24) describing the asymptotic behavior of f(r) and ^(r) 
can be used to determine the scattering matrix K directly from the quantum-defect matrix Y 

K = [Y(0)C„(0 + tt/2) - D„(0 + tt/2)]- 1 [Y(0)C,(0) - D,(0)] , (31) 

where O,(0) = {(SijO/fi)^)} and D v {4>) = {^jA^i)^)}- We note that in principle the 
parametrization in terms of the quantum-defect matrix Y(0) depends on the short-range phase 
that is fixed by the choice of the reference potentials. In contrast, the scattering matrix K will 
depend only on the actual scattering lengths associated with the scattering channels and should 
be independent of a particular set of reference functions that determine specific for that choice. 
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Figure 2. Energy-dependent scattering length a(k) as a function of the short-range phase <fi for 
different values of the collision energy. 



3.2. Scattering states and semiclassical approximation 

Using Eq.(26) and the small-g expansions presented in Appendix Appendix A. 4 one can derive the 
well known threshold behavior of the phase shifts in the polarization potential [23,30] 



Here, we have used the relation (23) between the short-range phase ip and the scattering 
length a of the reference potentials, in order to express all the formulas in terms of a. Because of 
the long-range character of the polarization potential, the scattering length can be defined only for 
I = 0, while for i > the threshold behavior of the phase shifts is dominated by the contribution of 
the polarization potential. Fig. 2 shows the energy-dependent scattering length a(k) = —tan^ /k 
as a function of for diferent values of the energy. The case k = corresponds to a = R* cot 4>. 

At sufficiently high energies the scattering can be described using the semiclassical 
approximation. To identify the crossover from the quantum to the semiclassical regime, we plot in 
Fig. 3 the quantity X'(r)/2-K, where A(r) is the local de Broglie wavelength. A necessary condition 
for the applicability of the WKB approximation, expressed in terms of the local wavelength, is 
A'(r) <C 2ir. We observe that this condition is first violated at distances comparable to R*, while 
the WKB approximation remains valid at small and large distances. We have verified that to a 
good approximation the wave function can be calculated within the semiclassical approximation 
when X'(r)/2iT < 1/2, a condition fulfilled for energies E > 25E*. Figs. 4 and 5 compare predictions 
of the WKB method with the exact wave functions f(r) for two values of energy: E = E* and 
E = 100E*. In addition they show the small-r and large-r asymptotic formulas (22a) and (24a), 
respectively, and derivative of the local wavelength: A'(r)/27r. 




3T + 0{q") 



t > 0, 



(32) 



(33) 
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Figure 3. Derivative of the local de Broglie wavelength A(r) giving the condition for the 
applicability of the WKB approximation plotted for different energies. 




Figure 4. Exact wave function (black solid line) compared with the WKB wave function (red 
solid line) and with small and large r asymptotic approximations, given by Eqs. (22a) and (24a), 
respectively (blue dotted and dashed lines), for energy kR* = 1 and = 0. The figure also shows 
the derivative of the local wavelength A(r) (green dot-dashed line), representing the condition for 
the applicability of the WKB approximation. 
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Figure 5. The same as Fig. 4 but for energy E = 100E* . 
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Figure 6. Energies of bound states as a function of the short-range phase </> for few lowest partial 
waves. 



3.3. Bound states 

For negative energies the exponentially decaying solution is given by a linear combination of f(r) 
and g{r) [see Eq. (9c)] with v given by (29). In the quantum-defect approach the bound state 
spectrum is determined by the condition Y + tan v(E) = 0, which is a single-channel version of 
Eq. (15). In the presence of a single channel one can take Y = 0, which is satisfied when the 
reference potential is chosen to reproduce the same scattering length as the real physical potential. 
This yields the condition tan v(E) = 0. In this case the bound state wave function <p{r) of Eq. (9c) 
contains only the f(r) component, with the short-range phase given in terms of the scattering 
length by Eq. (23). The condition determining the energy of a bound state can be rewritten as 

cos [§ v(E) + - f J V 16 / 

with S U (E) defined in the Appendix Appendix A. 3. Fig.6 shows energies of the bound states 
versus the short-range phase <fi for the lowest partial waves. We observe that for i even the bound 
states disappear at the threshold at <p = (a = ±oo), while for odd i this happens for <f> — tt/2 
(a = 0). Thus, for a polarization potential there is an I = 2 periodicity in the values of <fi at 
the threshold, which is reminiscent of I = 4 periodicity for van der Waals potential [33]. The 
energies of the bound states for a = ±oo and for different partial waves are shown in Fig. (7). 
For even I this defines characteristic energy bins, which determine the positions of the last bound 
states in the spectrum. For instance, the last s-wave bound state is located in the energy range 
E/E* = [-106,0] 

Applying small q expansions of v and S u (see Appendix A) one can expand the right-hand-side 
of Eq. (34), which yields for s- waves 

I = k--k 2 + C(k 3 ), £ = 0, (35) 
a 3 

and for p- waves 

V-^ + eV), £ = i, (36) 



15 
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Figure 8. Energies of s-wave bound states versus the inverse scattering length R* /a. Exact results 
of Eq. (34) (black solid line) are compared with the universal law E/E* = —(R* /a) 2 (red dotted 
line), and with the approximation (37) including terms up to the order of 1/a 3 (red dashed line). 



with k = y/— E. Inverting the former equation, we can write the energy of s-wave bound states in 
powers of 1/a 

E =-h +2 ih + °^> e=0 - (37) 

The first term on the right-hand-side is the universal energy of a weakly bound state, the higher 
order term represents the correction that is specific for 1/r 4 potential. Fig. 8 shows the binding 
energy for an s-wave bound state versus the inverse of the scattering length. We observe that the 
range of the applicability of the universal formula E = —1/a 2 is very narrow. Inclusion of the 
higher order correction in 1/a given by Eq. (37) improves the agreement with the exact one. In 
fact the approximation (37) works reasonably for \E\ < 0.5E*. 

3.4. Behavior of MQDT functions for polarization potential 

Figs. 9-11 show the C(E), ta,nX(E) and tanv(E) functions for angular momenta £ = and 2 
determined from Eqs. (27)-(29) for different values of the s-wave scattering length a i.e., choice 
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Figure 9. C(E) function versus the energy for partial waves I = (upper panel) and I = 2 
(bottom panel), calculated for different values of the scattering length a. 



of in Fig. 2. The calculations are performed under the assumption that the short range phase 
does not depend on £, such that all the MQDT functions can be parameterized with a single 
^-independent 0. Substituting the small-g expansions given in Appendix (Appendix A. 4) into 
Eqs. (27) and (28) one obtains 

C 2 (E) E ^° + ££-±ll_ sin 2 U + e *)(t\ 2 ' +1 + S ifi qcos 2 + 0(q 2 ), (38a) 



tan AGE) - cot U + £-) + 5 ifl q 2 ^- + 0(q 3 ), (38b) 



2 J ''"^ sin 2 iff 

.. \ . ft(fttan0 — 1) ttk 2 , n , ,. 

ta^W ~ ta, ^ + 1-) + 5 m > - s{( _ m + Otf. (38c) 

Here, 5e,o = 1 for £ = and is zero otherwise. At sufficiently large energy, the semiclassical solution 
begins to be valid at all distances, implying C(E) — > 1 and tanA(-E) — > 0. One can verify that the 
MQDT functions presented in Figs. 9-11 exhibit the small-g and large-g asymptotic behavior given 
by Eqs. (38). For nonzero angular momenta one can observe the existence of sharp peaks in the 
function C(E), corresponding to shape resonances, that appear due to the presence of quasi-bound 
states behind the centrifugal barrier. 

4. Frame transformation 



In the ultracold regime variations of the total energy E are typically much smaller than the depth 
of the potential at short range, where the matrix Y is defined. Therefore, it is justified to neglect 
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Figure 11. tanu{E) function versus the energy for partial waves I = (upper panel) and I = 
(bottom panel), calculated for different values of the scattering length a. 
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the dependence of ~Y(E) on energy and to set ~Y(E) = Y. In this way, the matrix Y can be 
determined at a single value of energy, for instance from the reaction or scattering matrices K 
or S by applying formulas (12) and (13), respectively. An alternative way is to apply the frame 
transformation technique [26,34] that provides a very convenient way of parameterizing the short 
range matrix Y in terms of few parameters, e.g. scattering lengths. 

In this paper we focus on the collisions of an alkali atom with an alkali-earth ion in their 
electronic ground states. Hence, the asymptotic channel states can be characterized by the 
hyperfme quantum numbers fi^m^ and f 2 ,rrif 2 for ion and atom, respectively, and by the angular- 
momentum quantum numbers i and of the relative motion of the atom and ion center of masses. 
We label the internal quantum numbers in the asymptotic basis by a = {fif-^Tn^m^}, hence the 
asymptotic channels are characterized by i = \a£rne). At short distances it is more convenient to 
characterize the channel states in terms of the total electron spin S = Si + s 2 and total nuclear 
spin I = ii + i 2 , where Si, s 2 are electron spin of ion and atom respectively, and ii, i 2 denote 
their nuclear spins, respectively. This basis is characterized by two additional quantum numbers, 
the total hyperfme angular momentum F — f j — f 2 — I — S and its projection Mf on the axis of 
quantization. We denote the internal quantum numbers in the molecular basis by = {ISFMp}, 
so the molecular channel states are characterized by j = \j3lme). The frame transformation is a 
unitary transformation between channels a and j3, that can be written as 

U aP = {hf 2 m fl m f2 \ISFM F ) = Y,{™ h m h \FM F ){hf 2 \IS). (39) 

F 

Here, for simplicity we have omitted the quantum numbers I and mi that are conserved if one 
neglects weak dipolar interactions. Symbol {mfimj 2 \FM F ) stands for the usual Clebsch-Gordan 
coefficient 

h h f 



m fl m f2 \FM F )=(-iy^ +M -V2FTT[ n JI , , (40) 

mf 2 M F 



while 



k si fx 

(fiMIS) =V(2fi + + 1)(2/ + 1)(25 + 1) { i 2 s 2 f 2 \, (41) 

I S F 

is the transformation between f\f 2 and IS coupling schemes, given in terms of the Wigner 9j 
symbol. 

The knowledge of the analytical solutions for the polarization potential suggest that the 
reference potentials can be chosen to contain only the diagonal long-range part of the interaction 
matrix: 



Vi(r) = FT + *TJ J ~ ^ C 42 ) 



KHijii + 1) C 4 
2/xr 2 

According to Eqs. (22a) and (22b) the wave functions in pure polarization potential are 
singular at r — > 0. In this case the standard boundary conditions F(r) — > (r — > 0) that is 
imposed on the physically meaningful solution, has to be replaced by the boundary condition 
(22a) with some short range phase fa. 
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4-1. Angular-momentum-insensitive quantum- defect matrix 

The exchange interaction that mixes asymptotic channel states a takes place usually in some 
range of distances Rq of the order of few tens of a®. At such distances the interaction potential 
between species is much larger than the hyperfine or Zeeman splittings in the presence of an 
external magnetic field. Therefore it is convenient to define the matrix Y at distances r > Rq, 
where the exchange interaction is negligible, and Wij(r) = SijC^/r 4 for Rq < r <C R*. This 
approximation relies on the fact that for r <C R* one can safely neglect the centrifugal potential, 
the asymptotic kinetic energy and hyperfine splittings. On the other hand, it also ignores the 
higher order dispersion terms Cg/r 6 , Cg/r 8 , • • • , which only give relatively small corrections for the 
specific atom-ion system considered in this paper [10]. Hence, for Rq < r <C R* the single-channel 
wave function is given by the linear combination of the two WKB solutions (22). 

Based on the arbitrariness of the reference potentials, we can choose values 4>% a ^ our 
convenience. For the following choice of short range phases 

<t>i = 0, (43) 
the quantum-defect matrix takes a particularly simple form in the molecular (IS) basis 

Yffi = 6pr[a 8W ]- 1 t (44) 

where S(/3) = 0, 1 denotes the total electron spin in the channel (3, and ao = a s and a\ = a t are 
the singlet and triplet s-wave scattering lengths, respectively. By using the unitary transformation 
(39) we find the quantum-defect matrix Y in the basis of the asymptotic channel states 

Y = UY (/5) U f . (45) 

In the presence of an external magnetic field B, the total transformation matrix is a product 
of U and an additional unitary matrix Z(JB) relating the bare (B = 0) and dressed (B 7^ 0) channel 
states \Wp(B)): 

\* a (B)) = Y,Z aa >{B) \m a ,{B = 0)) . (46) 

a' 

In this case the quantum-defect matrix reads 

Y = Z( J B)UY (/5) U t Z t (5). (47) 

Apart from the shifts of the internal hyperfine states of the atom and ion, the magnetic field in 
general also affects the motion of the ion as a charged particle. This leads the Landau quantization 
of the ion motion [35] and can result in several interesting scattering effects resembling the collisions 
in the quasi-lD confinements [36]. There are, however, two limits when the ion cyclotron motion 
can be neglected. First, when the characteristic size of the Landau orbit = with 
Q = eB/rrii and e denoting the ion charge, is much smaller than R*. Second, in the case when 
the ion is confined in the RF trap, which is typically much tighter than the size of the Landau 
states. This second case is however beyond the scope of the present analysis, and in the following 
we assume the former condition to apply. 
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4-2. Corrections to quantum- defect matrix for higher partial waves 

According to the arguments presented in section 4, at distances Rq <C R* one can expect the 
dependence of the quantum-defect matrix Y on £ to be negligeable. However, as we will show 
later with the example of a Na atom and a Ca + ion, this approximation works well only for the 
few lowest partial waves. Already for £ > 4 one starts to observe some deviations of the quantum- 
defect parameter from the exact numerical solutions for a realistic potential. In order to improve 
the accuracy of the quantum-defect model we introduce a correction of the short-range phase 
for nonzero relative angular momenta £. Such correction can be obtained within the semiclassical 
theory, since the modification of arises at distances where the real potential Wa{r) differs from 
Vi(r), a region well described in the semiclassical approximation. We start from the semiclassical 
formula for the radial wave function fi(r) calculated in potential Wa(r) 

£(r)=^(r)-V 2 sinA(r), (48) 



where ki(r) = ^2fi(E°° — Wa{R)/h is the local wavevector, /3j(r) = f + J^dxk^x) is the WKB 
phase, and tt is the classical turning point. Here, we neglect the contribution from the asymptotic 
kinetic energy, taking E = E°°. Eq. (48) can be rewritten in the following way 

Ur) = hir)- 1 / 2 sin [&(r) + A,(r)] , (49) 

where 

Ai(r) = A(r)- A(r), (50) 

h(r) = v/2/i(£°° - Vi(R)/h and A 

( r ) = i + Ir dxki(x). At sufficiently large distances, where 
the potential takes its asymptotic form (3), but still within the semiclassical regime r C fi*, we 
require that the phases of fi{r) and of the solution fi(r) calculated for Vi(r) be equal. At large r, 
Aj(r) — > Aj(oo) = const, and this happens already within the semiclassical domain r C i?*, since 
the potentials Wa(r) and Vi(r) have the same long-range asymptotics. The ^-dependent shift can 
be obtained by comparing the Aj(oo) computed for £ = and £ > 

5^ = ^(00)^-^(00)1^. (51) 

In this formula, AjL is calculated for the actual value of £, while Aj|^ =0 represents a shift that 
can be incorporated into the definition of the short-range phase 0. By treating the centrifugal 
barrier in the expressions for fcj(r) and ki(r) as a perturbation to the interaction potential, we find 
a particularly simple result containing partial-wave quantum number £ only in the prefactor of an 
integral 

m + i)h r°°dx ( 1 1 \ 

where U^r) = Wu(r) - E?° - ^ + l)/(2/ir 2 ). 
5. Radiative charge transfer 

In this section we develop a quantum mechanical description of the radiative charge exchange 
process A + B + — > A + + B + hv with emission of a photon of frequency v. Since the probability 
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of charge exchange in heteronuclear atom-ion collisions is small, the quantum transition rates 
can be described in the distorted wave Born approximation (DWBA) [37] treating the molecular 
dynamics exactly and the interaction with the radiation field as a perturbation. If weak dipolar 
interactions are neglected, the molecular dynamics conserves the quantum numbers £ and of the 
orbital angular momentum of the atomic fragments both in the initial and final molecular states. 
Radiative transitions are induced to first order by the electric dipole coupling, provided that the 
variation |A£| = 1. 

We ignore at first nuclear spins, such that the collision only involves two coupled molecular 
channels. After angular integrals have been performed analytically using standard angular 
momentum techniques, the radiative charge-exchange decay rate reads (see, e.g., Ref. [38]): 



This equation is expressed in terms of reduced dipole matrix elements for free-free (E'£'\d(r)\E£) 
transitions, in which the atom and the ion remain unbound, and free-bound (v'£'\d(r)\E£) 
transitions, in which the colliding pair forms a molecular ion. 

The \E£) scattering state is energy- normalized with incoming wave boundary conditions in 
the entrance potential of AB + state. The primed quantities \E'£') and \v'£') represent respectively 
energy-normalized scattering states with outgoing wave boundary conditions and unit normalized 
bound states of the A + B molecular complex. Energy conservation requires hv = E — E' + 5 and 
hv v i = E — E v > + 5 for free-free and free-bound transitions, where 5 is the difference of ion and 
atom ionization potentials. 

While this approach is physically transparent and has the advantage of providing detailed 
information on the products of the charge-exchange process, the generation of dipole moments 
as a function of photon energy and of free-bound matrix elements for all ro-vibrational levels 
is computationally costly. A significantly simpler approach discussed in [39] approximates the 
sum over all continuum and bound transitions with a simple average of a space-varying Einstein 
coefficient for spontaneous emission over the initial scattering wavefunction 



Eq. (54) is inspired by the exact quantum-mechanical closure relation over scattering and 
bound states in the final channel, and by the classical Frank-Condon principle relating the 
potential energy difference with the energy of a photon emitted at interatomic separation r, 
hv{r) = Wi(r) — Wf(r) [39]. In the next subsection we cast this equation in a particularly 
expressive form using a combination of semiclassical and MCQDT approaches. 



A(E) = ^{£/ d (H J [{£ + 1) \(E'£ + l|rf(r)|^)| 2 + £\(E'£ - l|rf(r)|^)| 2 ] 
+ [(£ + l)\(v'£ + l\d(r)\E£)\ 2 + £\(v'£-l\d(r)\E£)\ 2 ] }. (53) 




(54) 



with 




(55) 
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5.1. Semiclassical model 

The radiative charge transfer takes places at relatively short distances, where the electronic wave 
functions of the atom and ion start to overlap and the dipole matrix element is non negligible. At 
such distances one can calculate the matrix element of A(r) in the semiclassical approximation 

(E£\A{r)\E£) = lc~ 2 {E,£) [°° dr^f\. (56) 

Here, -R m ; n denotes the classical turning point in the potential U (r) of the entrance channel, and 




v(r) = \l-\lE-U(r)- h2{i 2 +J )2 (57) 

is the classical velocity of a particle in the entrance channel at the distance r. The MQDT function 
C~ 2 (E,£) provides the proper scaling of the semiclassical wave function at short distances, with 
respect to its long-range asymptotic behavior, given by the energy- normalized functions \E£). The 
total cross section for the charge exchange can be calculated from [40] 

a tr = ^hA(E). (58) 

In the regime of ultracold energies one can neglect the contribution of the energy E to v(r) at 
distances where the radiative charge transfer occurs. Similarly, one can also omit the contribution 
from the centrifugal barrier. In this approximation v(r) ~ y^—2U(r)/jj,, and the energy and 
angular momentum dependence enters only through the MQDT function C~ 2 (E,£) 

(E£\A(r)\E£)^±-C- 2 (E,£)P tr (59) 

Here, P tr = 2 dr A{r) /v(r) is the probability of the photon emission during a single collision, 
which in our approximation is a constant. This way the cross section for the charge transfer event 
o"tr can be written as 

a tr (E) = — P tl .^(2^+ 1)C^ 2 (P,£). (60) 

i 

It is instructive to investigate the classical, high energy limit of the charge transfer rate, where 
we expect the Langevin theory to be applicable. Starting from Eq. (54) with the semiclassical 
approximation (56), then setting C~ 2 (E,£) m 1 at E ^> E*, and replacing summation over £ by 
an integration over impact parameter b = {£ + we obtain [16] 

a g"*(E) = 2irJ2ii [°dbb [ dr . A ^ (61) 



J R 



mm 



E - U(r) - Eb 2 /r 2 



Eq. (61) can be further simplified, first by replacing the inner integral by a constant probability 
P t r 5 an d then by performing the remaining integral over b. We perform the integration over b to 
some maximal value of the impact parameter 6 max , given by the height of the centrifugal barrier 
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E max .(£) = \E*£ 2 (£ + l) 2 . Thus, the particle can penetrate the inner part of the potential only for 
E > E m&x , which imposes the upper bound on b: 6 max = ^J2R* /k. This yields for the cross section 

afr[E) = 2v^p tn (62) 
and for the charge transfer collision rate Ji tr = o~t r {E)v. 



K^ = ^P tI = 2nJ^P tI . (63) 

At large energies the charge transfer rate acquires a constant value. The latter formula agrees with 
predictions of the Langevin theory, assuming that all the classical trajectories that fall down on 
the scattering center lead to a reaction [22]. In our case the probability of a reaction in a single 
collision event is given by P tr , whereas in the case of homonuclear collisions, where the resonant 
charge transfer takes place P tT = \ [2]. 

It is convenient to express the charge transfer rate in the quantum regime in terms of K^ ass : 

K tT (E) = K* ass Q(E) (64) 
QW = 2jJrT,( 2i +l)C-*{E,£) (65) 

In this form the charge transfer rate is a product of a constant classical rate, that depends on 
the particular atom-ion system, and a quantum factor Q(E) that is universal and depends only 
on the specific atom-ion combination through the characteristic energy E* and the short-range 
phase 0. In this way the whole energy dependence will be given by Q(E). In the limit of high 
energy (E — > oo) Q(E) — > 1. In the limit of small energies (E — > 0), only s wave contributes and 
Q(E) |(1 + cot 2 0). 

In Fig. 12 we show Q(E) averaged over thermal distribution 

{Q{E))th ~ 0F(fc B T)3/ 2 J, dEQ(E)^e~ E ^ T (66) 

for selected values of the short range phase. We observe that even at large energies E ~ 10 A E*, 
(Q(E)) 7^ 1, due to the contribution of shape resonances. The detailed structure of resonances, 
however, is washed out by the thermal average. 



6. Results for 40 Ca + and 23 Na 

We now apply our quantum-defect model to 40 Ca + and 23 Na, a system whose ab initio potential 
energy curves are known [16]. This allows us to compare the predictions of our analytical approach 
with the full solution obtained by solving numerically the coupled-channel Schroedinger equation 
(1) with potentials of Ref. [16]. 
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Figure 12. Thermally averaged quantum contribution (Q(E)) t h to the charge transfer rate versus 
temperature for selected values of the short-range phase. 



Table 2. Scattering channels of 40 Ca + and 23 Na for Mp = 1/2 in the asymptotic representation. 
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6.1. Channel states 

The 40 Ca ion has vanishing nuclear spin (i\ — 0, s\ — 1/2). We label its ground-state 
sublevels in an external magnetic field in increasing energy \a\) = \f — 1/2, mf = —1/2) and 
\bi) = \f = 1/2, rrif = 1/2). The 23 Na has nuclear spin i 2 = 3/2, hence its hyperfine angular 
momentum can take the values / = 1 and 2. The hyperfine structure of 23 Na is shown in Fig. 13. 
The figure adopts the standard |a 2 ), \b 2 ), ■ ■ ■ notation to label states in increasing energy in the 
magnetic field, where in the weak field limit one has the identification \a2) = \f — ^,mf = 1), 
|6 2 ) = 1/ = 1, irif — 0), ■ • ■ . In the presence of an external magnetic field, only the projection Mj 
of the total angular momentum J = f x + f 2 + 1 is conserved during a collision. However, if we ignore 
small anisotropic spin-spin interactions giving rise to the long-range dipole-dipole force, the states 
with different angular momentum 1 will not be coupled, and both t and mi will be conserved. 
Hence, Mp will also be conserved and we can restrict our discussion to subspaces of constant Mp. 
In most of our calculations we will consider collisions within the subblock of Mp = 1/2, which 
contains the lowest energy channel state \a\a2). For Mp = 1/2 there are four possible scattering 
channels, listed in Table 2 together with their threshold energies at B — 0. Table 3 lists the 
channel states in the molecular basis. 

The quantum-defect matrix Y in the molecular {IS} basis is diagonal and for the assumed 
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Table 3. Scattering channels of 40 Ca + and 23 Na for M F = 1/2 in the molecular {IS} 
representation. 
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Figure 13. Hyperfine structure of 23 Na. Zeeman sublevels versus magnetic field B are shown. 
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(67) 



channel numbering it reads 

\ 

(atY 1 

{atY 1 

{atY 1 J 

Using the frame transformation (39) one can easily find Y in the asymptotic channel representation 

/ll y/2 V3 a/2 \ 

ia r 8a r 8a c 

1 

8a c 
1 1 

\ 8a c 4a c 8a c a t 4a c / 

where l/a c = l/a s — l/a t is the coupling parameter characterizing the strength of coupling between 
channels. We note that for similar triplet and singlet scattering lengths l/a c = 0, the channels are 
uncoupled, and no interchannel resonances occur. 



(k 8a c 

8a c 

8a c 
V2 



8a c 

1 1 Vq 

a t 4a c 8a c 

_\/6 1 3 

8a c at 8a c 

1 Vg 

4a r 8a r 



(68) 
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In the presence of a magnetic field, the energies of the hyperfine states of Na atom are given 
by the Breit-Rabi formula [41] 



Ef:]\B) 



E 



+ (-l) J 




(69) 



with 



x 



gfi B Bm f 
l + 2i 



y =(-i) 2(l+m ' Wb 



m. 



l + 2i 



where g is Lande factor, and /i# is the Bohr magneton. By including in the Hamiltonian additional 
small couplings of the magnetic field to the nuclear spin, B(B) = B(0) + /i B B(gjj + <?ii), the 
effective ^-factor is given by g = — and the levels acquire an additional B-dependent shift 
AEfr^j(B) = gifisBrnf. Here, g^ and g\ are Lande factors for the total orbital angular momentum 
j = L + s and the nuclear spin i, respectively. In Ca ion, which has no hyperfine structure, the 
levels shift according to the standard Zeeman formula 



= fiBBgjmf 



(70) 



The transformation matrix Z from the bare (B = 0) to the dressed channel states (B ^ 0) 
can be easily found from the transformation matrices Z^ k \ connecting bare and dressed states of 
ion (k = 1) and atom (k — 2). For the ion we have simply = 1, while for the atom 
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(71) 



where f 2 = i 2 + 1/2 = 2, and w = y/y 2 + (x + Ej 2mn 
for the assumed channel numbering reads 



2 . Then the total transformation matrix Z 
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6.2. Bound states 

In this section we discuss properties of bound states in a single scattering channel. Fig. 14 shows a 
sample bound state spectrum for a 3 E + molecular potential. The numerical calculations based on 
the ab initio Na-Ca + potentials of Ref. [16] with the quantum-defect theory with and without the 
^-dependent correction (51) to the short-range phase are compared. For simplicity in the numerical 
calculation we neglect the higher-order dispersion terms in the potential, so to isolate the effect of 
the centrifugal barrier, in particular for high order partial waves. We note that the quantum-defect 
model assuming the same short-range phase for all partial waves starts to deviate already for 
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<j)ln 

Figure 14. Energies of the bound states versus the short-range phase 4> in the triplet channel a 3 S + 
for few lowest partial waves. The numerical results obtained by solving the Schrodingcr equation 
with ab-initio potential of Rcf. [16] (solid lines), are compared with predictions of MQDT assuming 
angular- momentum insensitive short-range phase (dotted lines), and with MQDT including the 
^-dependent correction (51) to the short-range phase (dashed lines). The inset zooms in the region 
around = 0. 

I = 3, whereas the inclusion of the correction (51) greatly improves the agreement between MQDT 
model and numerical solution. The inset shows magnification of the region around = 0. For 
I = 4 still some small discrepancy between the numerical and corrected MQDT results can be 
observed. Tha latter may originate from effects going beyond the WKB approximation, which was 
assumed to hold in the derivation of Eq. (51). 

6.3. Radiative charge transfer 

In the case of Na-Ca + collisions the charge transfer process occurs due to transitions from a 
singlet A y Tj channel to the absolute ground molecular term NaCa + (A 1 S) — > Na + Ca(X 1 E). For its 
description we first use the DWBA leading to Eq. (53), which allows the final states of the charge 
exchange process to be individually identified. In order to calculate the free-free transition dipole 
elements a two coupled channel model is setup, comprising Wi(r) = V [A 1 !]) and Wf(r) = V (AT X E) 
diagonal molecular potentials with I and l' centrifugal barriers, coupled by a £<i(r) term, with £ a 
parameter to be optimized. In a field-dressed approach, a vertical energy shift is introduced so to 
guarantee total, photons plus atoms, energy conservation. If the £ parameter is chosen so to insure 
validity of first order perturbation theory, in the DWBA the required reduced matrix element is 
simply proportional to the transition scattering matrix element s^ = s(E'£' 4— E£; £) obtained 
numerically 

s ? = -2iti£(E'l'\d{r)\E£). (73) 

The free-bound reduced dipole elements are obtained following a different method. As a first 
step, the vibrational wavefunctions of the X l H potential are generated for each angular momentum 
£' using a grid method. To this aim use of the step-adaptive approach of Ref. [42] is essential. 
In fact, due to the long-range character of the r -4 interaction, numerical convergence of near 
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Figure 15. The closest to threshold £ = vibrational wavefunction for s-wave scattering length 
a = 5R* (lower panel) showing the points of the adaptive numerical grid. Convergence of the 
molecular energy is also shown (upper panel) as a function of the numerical grid size for two values 
of the radial step Sr near the potential minimum (see text). 



threshold levels is attained for a grid extending to large distances. To keep the number of grid 
points reasonable (~ 5000), the local step size is made to vary over five orders of magnitude from 
the potential well to the asymptotic region using an appropriate scaling function [42]. Fig. 15 
shows for instance the last rotationless vibrational wave function for a = 5R*, as well as the 
convergence rate of the corresponding eigenenergy for two different step sizes 5r of the numerical 
grid near the potential minimum. Note that for this specific molecular level, convergence begins 
to be attained for a box of 10 5 ao size. The relative energy accuracy is on the order of 10~ 2 
for Sr = 0.2a , and increases by one order of magnitude for Sr = O.lao- As a second step, the 
overlap (v ; f\d(r)\E£) with the initial scattering wavefunction of angular momentum £ is computed 
using a standard propagation code. The number of needed £ values is determined by numerical 
convergence. 

In the alternative approach based on the approximate sum rule Eq. (54), the expectation value 
(E£\A(r)\E£) are extracted numerically from the elastic scattering matrix element = s(E£ 
Ei\ £) for a single channel potential Wi with angular momentum I perturbed by a ^Air) term. The 
scaling parameter £ has to be chosen so to be in the linear regime. With this proviso, letting sq be 
the scattering matrix for the unperturbed Wi potential, the needed matrix element is determined 
from the DWBA as 

sSfo - s o) = -2ititi(El\A{r)\El). (74) 

The prefactor Sq arises from the fact that the expectation value on the right hand side has to be 
evaluated between scattering states with outgoing boundary conditions whereas transition matrix 
elements in the DWBA present incoming wave boundary conditions in the exit channel. 

The charge transfer rates we obtain are shown in Fig. 16 as a function of the collision energy for 
a sample value of the singlet scattering length, a s = R* . The figure compares the numerical result 
calculated by summing contributions from all free-free and free-bound transitions Eq. (53) with 
the MQDT model assuming the semiclassical charge-transfer probability Eq. (64). In the former 
case we additionally plot separately contributions from free-free and free-bound transitions. The 
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Figure 16. Rates of the radiative charge transfer due to the transitions between A 1 E + and 
X 1 Tj + states versus energy for the single scattering length a s = R* . The full numerical calculation 
based on the Fermi golden rule Eq. (53) (triangles) is compared to the quantum-defect model 
assuming semiclassical description of the charge transfer process (solid line). In addition we show 
the contribution from free-free (squares) and free-bound (circles) transitions obtained from the 
numerical calculation, the thermal average obtained from the quantum-defect model (dashed line) 
and the charge transfer rate given by the Langevin theory (gray solid line) . 



MQDT calculation includes the ^-dependent correction (51). We observe that at energies larger 
than 10/iK, the rate of charge exchange exhibits several peaks due to the shape resonances. The 
MQDT model agrees well with the full numerical calculations, except in the mK regime where 
high order partial waves are contributing. In this range of energies the discrepancy is due to the 
corrections to the short-range phase for large I that are beyond the applicability of the semiclassical 
formula (51). Finally, the approximation based on Eq. (54) is in full agreement with the numerically 
exact result and is not shown in the figure. 

6.4- Population of vibrational states in the charge transfer process 

According to Fig. 16 the main outcome of the charge transfer process are molecular ions. It is 
interesting to analyze the dependence of the charge exchange rate on the vibrational quantum 
number of the molecule, i. e. the vibrational distribution of the product molecular ions. In the full 
quantum treatment based on the numerical calculation of initial and final wave functions, such 
distribution is proportional to the transition matrix elements in Eq. (53). In the MQDT approach 
it can be obtained using the semiclassical approximation, in analogy with the derivation presented 
in Section 5.1. We start from Eqs. (54) and (56), and obtain the probability of the charge transfer 
process normalized per vibrational quantum number in the final state 

^mkm^V' J+1 ^ E,J) - (75) 

Here, p(E) = 4§ is the density of states in the exit channel and V(r) = Wj(r) — W{(r) is the 
difference between interaction potentials of the entrance (initial) and the exit (final) channels. The 
distance r at which all r-dependent quantities are calculated is related to the final energy by the 
Frank-Condon principle: E + 5 — E' = hu(r) = V(r), with 5 = E?° — Ef° denoting the difference 



Multichannel quantum- defect theory for ultracold atom-ion collisions 



27 




4x10" 



-0.015 



-0.010 



-0.005 



0.000 



Figure 17. Rates of the radiative charge transfer due to the transitions between A 1 Y I + and 
X 1 ^ states versus energy of the vibrational state in the exit channel, for two values of the singlet 
scattering length a s (lines with symbols). The solid line shows the semiclassical result obtained 
from Eq. (75), while the dashed line is the same result multiplied by two, which gives approximately 
the amplitude of the quantum oscillations. 



of the dissociation energies of the entrance and exit channels, respectively. The density of states 
p(E) can be calculated from the LeRoy-Bernstein formula [43]. In the the case of r -4 potential it 
yields 



In Fig. (17) we show the charge transfer rate from A 1 S + channel to the bound states of the 
X 1 E + channel, calculated numerically in the DWBA. The numerical points are obtained for two 
values of the singlet scattering length: a s = ±R*. The solid line shows prediction of Eqs. (75), with 
the total cross section related to the probability A by Eq. (58). We note that the full quantum 
result predicts oscillations of the distribution over final states, whereas the semiclassical theory 
leads to a smooth behavior, which can be interpreted as the distribution averaged over the quantum 
oscillations. The origin of these oscillations in the quantum result is due to the effects of the phase 
matching between the entrance and the exit channel wave functions, which is not present in the 
semiclassical result assuming Frank-Condon approximation. 

6.5. Magnetic Feshbach resonances at zero energy 

In this section we analyze magnetic Feshbach resonances in the limit of zero energy when only 
s-wave scattering is present. We neglect the effects of the magnetic field on the translational ion 
motion. For the Na-Ca + system of interest, the radius of the Landau orbit ol becomes equal to 
R* at Bl = 1086(7. Hence, the presented analysis is valid as long as B <C B^. 

Figure 18 shows the variation of the s-wave scattering length versus magnetic field for two 
sample values of the singlet and triplet scattering lengths: a s = R* = 2081ao and a t = —R* = 
— 2081a , respectively. The resonances can be assigned to the particular bound state crossing the 
threshold of the a = \a1a2) channel, as one can verify by inspection of the bound state energies in 




(76) 
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the bottom panel. The resonances are labeled by capital letters from A to E. The middle panel 
shows a close-up of the bound-state spectrum just below threshold. Because of the long-range 
character of the polarization potential, the last s-wave bound state is always located relatively 
close to the threshold, with the binding energy Ef, < 106E* = 3.02 MHz. As one can observe in 
the figure the last bound state in the open channel can be strongly coupled to the other bound 
states crossing the threshold, giving rise to relatively strong avoided crossing, as in the case of 
resonance E. In contrast, very close to threshold the energy of the bound state bends and follows 
the universal behavior E = —h 2 (2fia 2 ). 

An approximate but highly accurate description of the Feshbach resonances can be developed 
using a two channel description based on CI (configuration interaction) model of a Feshbach 
resonance, where a single closed channel represents the effects of all the closed channels contributing 
to a resonance [44]. In the case of a single open and a single closed channel, the application of 
MQDT is straightforward, and after some simple algebra (see for instance [45]), one obtain the 
following expression for the phase shift £(-E, I) in the open channel [46] 

£(EJ) J) -tan- 1 | 2 ° } E ^\ ) . (77) 

The first term £bg(-^ ; is the background phase shift, describing the scattering from the open 
channel only and incorporating effectively the influence of the closed channels for magnetic fields 
far from the resonance. The second term describes the resonant contribution resulting from a 
bound state in the closed channel with energy E n , crossing the threshold of the open channel. 
An energy- dependent width of the resonance is given by a constant width T multiplied by the 
MQDT function C~ 2 (E,£), which accounts for a proper threshold behavior as k — > 0. In the 
two-channel model one assumes that the energy of a bound state varies approximately linearly 
with the magnetic field B: 

E n (B)=8fi(B-B n ), (78) 

where B n is the magnetic field at which the bound state crosses the threshold of the open channel, 
and 5fi is the difference of magnetic moments between the open and closed channel. The second 
MQDT function tanX(E) in the denominator of £ res , describes the shift of the resonance position 
from the bare value B n , which is due to the coupling between the open and closed channel. 

In the following we will focus on s-wave Feshbach resonances, considering the zero-energy 
limit. According to the results of Sections 3.2 and 3.4, for I = the phase shift and MQDT 
functions exhibit the following threshold behavior 

tan£(£) -ha - \{kR*) 2 , (79) 

C~ 2 (E) E ~° kR*(l + s 2 ) - (fciT)V(l + s 2 ), (80) 
tanA(E) E ~° -s + (kR*) 2 s(l + s 2 ), (81) 

where s = a/R*. The connection to the standard theory of magnetic Feshbach resonances can be 
done by introducing the width of magnetic Feshbach resonance A: 

l im £ C 2 ^ E \ = Sfi A, (82) 
2tan£ bg (£) r K J 
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and the resonance position B , that is shifted from B n due to the coupling between the open and 
closed channel 

r 

B = B n + — lim tan X(E). (83) 

Now, substituting Eqs. (78), (82), (83) into (77) in the limit of zero kinetic energy one retrieves 
the standard expression 



a(B) 



A 

1 - 



4) 



B-B 

Making use of the MQDT expansions Eqs. (79)-(81) one can express T and B n in terms of 
parameters a^ g , A, B , S/j, which can be directly measured in experiments 

B n = B + -^VA, T = t^H//A. (85) 

Here, s hg = a hg /R*. 

Another important parameter characterizing Feshbach resonances is the fraction of the closed 
channel in the weakly bound molecular state at large and positive values of the scattering length 
[47] 

Ba-B 



Z(B) - 1 ^ - 1 



A 



6) 



Sfj, dB ( 

The parameter ( describes the range of magnetic fields expressed as a fraction of the resonance 
width A over which the resonance exhibit the universal properties and the occupation of the 
closed channel remains small. Sufficiently close to the resonance, the binding energy is given by 
the universal formula E = —h 2 (2fia 2 ), which leads to [47] 



•s 



C=f-f?- (87) 

For entrance channel dominated resonance £ ^> 1 and Z(B) remain small for detuning of the order 
of A. In the opposite case (<1, the Feshbach resonance is called closed-channel dominated, and 
the universal regime where the energy-independent formula (84) is applicable is very narrow. 

We have fitted the universal formula (84) to Feshbach resonances presented in Fig. 18. The 
results are summarized in Table. 4. The difference of the magnetic moments has been determined 
from the bound state spectrum (lower panel of Fig. 18) by extracting the linear slope for the 
molecular states giving rise to the resonances. Finally, the parameter £ has been calculated from 
Eq. (87). We note that only the two first Feshbach resonances, occurring at relatively small 
magnetic fields, are non universal, while the remaining ones are relatively broad and entrance- 
channel dominated. 



\6»A\ 



6. 6. Magnetic Feshbach resonances at finite energies 

We analyze here the dependence of elastic and inelastic rates on the magnetic field at finite collision 
energy, where several partial waves play in general an important role in the collision physics. To this 
aim we use both our MQDT model and numerical calculations based on the close-coupled (CC) 
Schroedinger equation. In order to calculate the charge exchange rates from the CC solutions 
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Figure 18. Scattering length (upper panel) and energies of the s-wave bound states (middle 
and lower panel) versus the magnetic field strength B for collisions of Na with 40 Ca + calculated 
using MQDT model for some typical singlet and triplet scattering lengths: a s = R* = 2081ao 
and a t = —R* = — 2081ao, respectively. The capital letters in the upper panel label the Feshbach 
resonances. 



Table 4. Parameters of Feshbach resonances shown in Fig. 18: resonance position Bq, resonance 
width A, background scattering length etb g , difference of magnetic moments between the open and 
closed channel o/x, and parameter £ characterizing the fraction of the closed channel (see Eq. (87)). 
The labels in the first column enumerate the resonances in accordance with markings in Fig. 18. 
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0.322 


-0.000417 


-2019 
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B 
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4.46 
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-1787 


3.15 


79.6 
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we have generalized the approximate closure relation Eq. (54) to the multichannel case, including 
hyperfine effects and the external magnetic field. The charge-exchange process in the MQDT model 
has been described using the semiclassical approximation Eq. (60) where the MQDT function 
C~ 2 (E) is replaced by the amplitude of the wave function in the singlet A l T> channel. Details are 
provided in Appendix Appendix B. 

Figs. 19, 20 and 21 show the elastic and charge transfer rates at three collision energies: 
E = InK, E = 1/iK and E = lmK for some typical values of the singlet and triplet scattering 
length, a s = R* = 2081ao and a t = —R*. Contributions from the lowest 2, 9 and 21 partial waves 
are included respectively at these three energies. At InK resonances appear mostly in the s-wave 
and are relatively broad. In addition one can observe few narrow resonances in the p-wave channel 
occurring at B = 27, 84.7 and 182G. In this range of temperatures the MQDT model is extremely 
accurate and agrees perfectly with the full numerical CC calculations. 

At higher energy E = 1/iK the agreement of the analytic model with the numerical solution is 
still very good. At this energy the resonance peaks arise from the 4 lowest partial waves. The arrow 
in the bottom panel of Fig. 20 indicates the I = 3 Feshbach resonance. It is the only resonance at 
E = 1/iK whose position is not well predicted by the analytical MQDT model. We note that the 
charge transfer rates exhibit more resonance peaks than the elastic rates. This general behavior can 
be qualitatively understood by analyzing the number of partial waves contributing to the elastic 
and inelastic processes. The elastic rates are dominated by the reflection from the long-range 
polarization potential and their contribution to the cross section decay as (2£+ 1) sin 2 ^(/c) ~ l/f 5 
at large I (see Eq. (33)). In contrast the charge exchange process must involve the tunneling 
through the centrifugal barrier, and its probability decays exponentially with I. Hence the number 
of partial waves contributing to the inelastic process is much smaller than for elastic scattering, 
and the inelastic rates are more sensitive to scattering resonances. On the other hand, the narrow 
resonances from high-order partial waves that appear in the elastic rates are less pronounced due 
to the strong background arising from reflection on the long-range potential. 

Finally, at the highest energy considered E = lmK the resonances are narrower and have 
smaller amplitude. This effect arises due to the large number of partial waves contributing to the 
scattering, which have tendency to wash out the resonance structures. We observe that at lmK 
the MQDT model basically follows the magnetic field dependence of the exact numerical rates, 
but it predicts accurately only the resonances associated to the lowest partial waves. We note 
that inclusion of the thermal averaging washes out the resonance structure for the elastic rates, 
while the charge-transfer rates still exhibit some resonance peaks. This is again due to the vastly 
different numbers of partial waves contributing to the elastic and inelastic collision processes. 

7. Conclusions 

Summarizing, we have developed a quantum-defect model for ultracold atom-ion collisions. The 
model was applied to the reference system composed by a 40 Ca + ion and a 23 Na atom, and 
its predictions were thoroughly verified by comparison with numerical close-coupled calculations 
using ah initio potential energy curves. Our model is based on the multichannel quantum-defect 
formalism, where the quantum-defect parameters are defined in terms of the analytic solutions 
for r -4 polarization potentials. Use of a frame-transformation allows us to reduce the number of 
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Figure 19. Elastic (upper panel) and charge exchange (lower panel) collision rates for collisions of 
Na and 4 0Ca + versus magnetic field calculated at energy E = InK for singlet and triplet scattering 
lengths a s = R* = 2081ao and at = —R* = — 2081ao. The figure compares the CC numerical 
calculations (black solid) and the quantum-defect model (red dashed). 
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Figure 20. Same as Fig. 19 but for collision energy E = The arrow in the lower panel 

indicates the position of the I = 3 Feshbach resonance (see text for details). 
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Figure 21. Same as Fig. 19 but for collision energy E = lmK. The thick blue line represents the 
MQDT result averaged over a thermal distribution with temperature of lmK. 



short-range parameters to essentially singlet and triplet scattering lengths only. Since for atom- 
ion systems of experimental interest the values of the singlet a s and triplet at scattering lengths 
are not yet known, in our calculations we have assumed typical scattering lengths of the order 
of the characteristic length R*. Once a s and a t will be measured our model could be readily 
applied to obtain all the basic collisional properties in the ultracold domain, including accurate 
positions of Feshbach and shape resonances. Application of our theory to other atom-ion systems 
is straightforward, amounting to a simple change of the scattering lengths and to the use of new 
characteristic parameters R* and E*, which are determined by the atomic polarizability and the 
reduced mass. 

In our studies we considered only two-body collisions in free space ignoring possible effects of 
the trapping potential, which is present in realistic systems. This can be of particular importance 
for the ions, that acquire in the presence of a time- dependent radio-frequency potential a small 
amplitude high frequency motion known as micromotion. This effect is important for current 
experiments, leading on one side to a significant loss of atoms, on the other potentially preventing 
sympathetic cooling of the ions to the ground-state of the trapping potential [48]. Moreover, 
the presence of a tight ion trap with characteristic size smaller than R* modifies the long-range 
asymptotics of the atom-ion wave functions and results in principle in the coupling of relative and 
center-of-mass motion. In addition, in the presence of a magnetic field the charged ions describe 
cyclotron orbits, an additional confinement effect which can lead to the appearance of scattering 
resonances [36] and may affect the actual positions of Feshbach resonances. We are currently 
investigating these issues. 
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Appendix A. Properties of analytical solutions for r potential 

Appendix A.l. Basic derivation 

We start from the Mathieu's equation of the imaginary argument (18) 

d 2 ^ 



dz 2 



[a - 2gcosh2;z] ij) = 0. (A.l) 



where a = (I + |) 2 and q = \JE. It is convenient to look for the solution of (18) in the following 
form 

oo 

M v {z) = c n {v)e^\ (A.2) 

n=— oo 

where v is the characteristic exponent. Substituting the ansatz (A.2) into (18) we obtain the 
following recurrence relation 

[(2n + u) 2 -a]c n + g(c n _! + c n+l ) = 0. (A.3) 

The three-term recurrence can be solved in terms of the continued fractions. In analogy to the 
solution of r -6 potentials [49] we subsitute 



q y r(^ + i)r(^ + i 



(A.4) 



r(^-n)r(^- B 



b~ (A.5) 



for n > 0. In the case of n = we have c = b$ = b$ . Now the recurrence relation (A.3) can be 
written as 

h+ -h+ = q n+l CA 6) 

n [(2n + 2 + v) 2 - a] \{2n + v) 2 - a] ' 1 

n 2 b~ 

u- _ h- _ q Vj ( A J\ 

n [(2n + 2 - z/) 2 - a] [(2n - z/) 2 - a] ' 1 ' J 

Finally we substitute fa+ = b^/b^^ and /i~ = b~/b~_ x , which yields the continued fractions 
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To find values of the coefficients c n it is sufficient to set /i+ = 1 and h m = 1 for some, sufficiently 
large m and calculate /i+ and h~ up to n = 1 using (A. 8) and (A. 9). Then 

bt = h+h+_ 1 ...h+c , (A.10) 

K = h n h n-l--- h I C 0, (A. 11) 

and coefficients c n can be obtained from Eqs. (A.4)-(A.5). Characteristic exponent v has to 
determined from Eq. (A. 3) with n = 0: 

h±(y) K{v) 



2 2 

z/ — a — q 



0. 



(z/ + 2) 2 -a (z/-2) 2 - a/ 
In numerical calculations it is more convenient to find v from equation [32] 

COS TXV = 1 — A(l — COS7T\/a) 

where A is an infinite determinant (independent of v) 
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(A.12) 



(A.13) 



(A.14) 



with 7 n = g/(4n 2 — a). Typically determinant A converges very fast and to calculate A it is enough 
to take relatively small matrices. 



Appendix A. 2. Asymptotic expansions for large arguments 

To derive asymptotic expansion of M u (z) for z^oowe observe that the leading contribution to 
the sum (A. 2) comes from the terms with large n. We neglect contribution of terms with n < 
and apply the following approximation for the terms with n > 



q\ n 



+ 1 



4 

where fc+ = lim^oo 6+. This yields 

' v - 



v—ifa 



n r 



+ 1 + n 



(A.15) 



M v {z) 



b + r 



i r 



V 



(A.16) 



where J v denotes Bessel function. Now the asymptotic behavior for large z can be easily obtained 
from the well-known asymptotic expansions of the Bessel funtions. 
In similar way we obtain behavior for large and negative z 



r i 



V 



^] J-u(y/ge- 



(A.17) 



where = lim n ^ 0O b n . 
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Appendix A. 3. Two linearly independent solutions 

As the two linearly independent solutions of (18) we can take M v (z) and M_ u (z) f. It is convenient 
for our purposes to define the following two linearly independent solutions in initial variable r 

T u {r) = w[v)M v [In (y/qr)] y^F, (A.18) 

T- V (r) = w(-v)M_ u [In {Jqr)\ (A.19) 

where 

w(u) — \ ; : „ ' \, 7^ TTV A.20 

Using Eqs. (A. 16) and (A. 17) one can easily work out asymptotic behavior of T v (r). For small 
r we obtain 



T v ( r )^r C os[^- + -v--), (A.21) 



1 7T 7T 



whereas for large r we get 



where 

( )_ feMr(^ + f)r(l-^) (A ' 23) 

Appendix A. 4- Expansions for small q 

The small-g expansion of the characteristic exponent v can be obtained from Eq. (A. 12) 

-'4- 4( i -l)( i + |)( i + l) + °^- < A ' 24 > 

Utilizing this results and applying (A.8)-(A.9) and (A.lO)-(A.ll) we calculate expansions of h^iv) 
and 6~ [y) 

^)-^ (,+ ^g- + s:? +i) v + ^). <-) 

where ^(cc) denotes the digamma function. 

Finally expansion of S[y) can be obtained from definition (A. 23) where we substitute (A. 25), 
(A. 26) and expand the Gamma functions. This yields 

( ' ~ IXi - L 1 4(/-i) 2 (/ +!)((+ |) 2 «+°<«>]- (A.27) 

| It is easy to oberve from Eq. (A. 3) that if v is a characteristic exponent then — v must be also a characteristic 
exponent and c n (— i/) = c_ n {v) 



Multichannel quantum- defect theory for ultracold atom-ion collisions 



37 



Appendix B. Multichannel calculations of the radiative charge transfer 

Equation (54) based on the approximated closure relation can be easily generalized it to the 
multichannel case including hyperfine effects and an external magnetic field. In this case the 
initial state for atoms incoming in the dressed channel state a of Eq. (46) will be labeled \E£a). 
The reduced dipole moment becomes a matrix, with only nonvanishing diagonal elements if the 
total electron and nuclear spin {SI} representation is used. 

In the presence of a magnetic field close-coupled scattering equations are numerically solved 
with a matrix perturbation £A(r). The emission rate for atoms incoming at collision energy E is 
expressed in the DWBA in terms of the open-open scattering matrix elements S£p a = s(E£/3 <— 
E£a;£) with and without field perturbation 

- so/3a) = -2m^(E£a\A(r)\E£a). (B.l) 

P 

As in the single-channel case, right multiplication by the Sq matrix enforces the correct boundary 
conditions, and symmetry of the scattering matrix {s^ a p = S£p a ) resulting from time-reversal 
invariance has been used. 

Calculation of the charge-exchange rates in the MQDT approach can be done using the 
semiclassical formula Eq. (60), with the function C~ 1 (E,£) replaced by the amplitude A{ IS \e,£) 
of the singlet component of the multichannel wave function at short range 

°tr(E) = ^P tT J2(2£ + l)\A[ IS \E,£)\' 2 (B.2) 
i 

(IS) 

In order to calculate A s we first analyze the open-open block of the multichannel wave function 
at large distances 

foo(r) + g o(r)Y J A OQ (B.3) 



r) — > 



Using relations between short-range and long-range normalized solutions, Eq. (B.3) can be 
rewritten as 



/ s r— >oo 

r — ► 



foo(r) + goo (r)R(£)] C(E) [1 - tan X(E) Y 00 ] A 00 (B.4) 



with R(-E') given by formula (11) applied for the open-open block of the renormalized quantum- 
defect matrix Y OQ . The constant matrix A OQ is fixed by the boundary conditions at r — > oo. With 
the following choice of A QO , 

A 00 = [1 - tan X(E)Y 00 ] ^ C(E)' 1 [1 - iR(£')] _1 e* (B.5) 

the wave function has a normalization corresponding to a unit flux of incoming particles 

Foo(r) ™ - [R^ 2 \kr) + H«(fcr)S] . (B.6) 
2 



Here, H\-\kr) — > dijie l ^ kiT ti7V ^/y/ki and H^'ikr) = H>j (kr) are functions exhibiting 

asymptotic behavior associated with the spherical Hankel functions hf\kr) and hg(kr), 
respectively. 
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The total wave function expressed in terms of the short-range normalized solutions reads 



F(r) 



with 



A 



f (r) + g(r)Y 



A 



- [Y 00 + tan Y co A c 



(B.7) 



(B.8) 



chosen in such a way that closed channel wave function at large distance is proportional to the 
exponentially decaying solution <f)i(r), Eq. (9c). Applying the frame transformation yields the 
multichannel amplitude of the wave function in the molecular basis 



A (7S) = (Z(i?)Ucosx) _1 A 



(B.9) 



where xpp' = fipp ,( fis(p) is the diagonal matrix containing short-range phases of the singlet (fio and 
the triplet <f>\ potentials, with S(/3) denoting the total electron spin in the channel (3. 
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